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Abstract 

A novel optimization method is proposed to minimize a convex function subject to bilinear 
matrix inequality (BMI) constraints. The key idea is to decompose the bilinear mapping as a 
iy) ■ difference between two positive semidefinite convex mappings. At each iteration of the algorithm 

the concave part is linearized, leading to a convex subproblem.Applications to various output 
feedback controller synthesis problems are presented. In these applications the subproblem 
in each iteration step can be turned into a convex optimization problem with linear matrix 
inequality (LMI) constraints. The performance of the algorithm has been benchmarked on the 



O . data from COMPl c ib library. 



(-H t Keywords: Static feedback controller design, linear time-invariant system, bilinear matrix 

inequality, semidefinite programming, convex-concave decomposition. 



1 Introduction 

Optimization involving matrix constraints have broad interest and applications in static state/output 
feedback controller design, robust stability of systems, topology optimization (see, e.g. [51 151 I2"T1 [T§] ) . 
Many problems in these fields can be reformulated as an optimization problem of linear matrix 
inequality (LMI) constraints [5j [21] which can be solved efficiently and reliably by means of interior 
point methods for semidefinite programming (SDP) [31 [55] and efficient open-source software tools 
0^ ' such as Sedumi [31] . SDPT3 [33 . However, solving optimization problems involving nonlinear matrix 

inequality constraints is still a big challenge in practice. The methods and algorithms for nonlinear 
matrix constrained optimization problems are still limited [9] 1 1 1 1 119) . 

In control theory, many problems related to the design of a reduced-order controller can be 
conveniently reformulated as a feasibility problem or an optimization problem with bilinear matrix 
inequality (BMI) constraints by means of, for instance, Lyapunov's theory. The BMI constraints 
make the problems much more difficult than the LMI ones due to their nonconvexity and possible 
nonsmoothness. It has been shown in [3] that the optimization problems involving BMI are NP-hard. 
Several approaches to solve optimization problems with BMI constraints have been proposed. For 
instance, Goh et al [12] considered problems in robust control by means of BMI optimization using 
global optimization methods. Hoi et al in [16] proposed to used a sum-of-squares approach to fixed 
order TJ-infinity synthesis. Apkarian and Tuan [5] proposed local and global methods for solving 
BMIs also based on techniques of global optimization. These authors further considered this problem 
by proposing parametric formulations and difference of two convex functions (DC) programming 
approaches. A similar approach can be found in pQ. However, finding a global optimum is in 
general impractical while global optimization methods are usually recommended to a low dimensional 



* Department of Electrical Engineering (ESAT-SCD) and Optimization in Engineering Center (OPTEC), K.U. 
Leuven, Kasteelpark Arenberg 10, B-3001 Leuven, Belgium. Email: {quoc. trandinh, moritz. diehl}® esat.kuleuven.be 

^Department of Computer Science. Katholieke Universiteit Leuven, Leuven, Belgium. Email: {suat.gumussoy, 
wim. michiels}@cs. kuleuven.be 



1 



problem. Our method developed in this paper is classified as a local optimization method which 
aims as finding a local optimum based on solving a sequence of convex semidefinite programming 
problems. Sequential semidefinite programming method for nonlinear SDP and its application to 
robust control was considered by Fares et al in [10] . Thevenet et al [34] studied spectral SDP methods 
for solving problems involving BMI arising in controller design. Another approach is based on the 
fact that problems with BMI constraints can be reformulated as problems with LMI constraints with 
additional rank constraints. In [53] Orsi et al developed a Newton-like method for solving problems 
of this type. 

In this paper, we are interested in optimization problems arising in static output feedback con- 
troller design for a linear, time-invariant system of the form: 



where x £ W 1 is state vector, w £ 1"" is the performance input, u £ W n " u is input vector, z £ l" z is 
the performance output, y £ W ly is physical output vector, A £ M. nxn is state matrix, B £ M nxn * 
is input matrix and C £ R n » x ™ is the output matrix. Using a static feedback controller of the form 
u = Fy with F £ M. nuXrly , we can write the closed- loop system as follows: 



The stabilization, , "Hoo optimization and other control problems for this closed- loop system will 
be considered. 

Contribution. Many control problems can be expressed as optimization problems of BMI con- 
straints and these optimization problems can conveniently be reformulated as optimization problems 
of difference of two positive semidefinite convex (psd-convex) mappings (or convex-concave decom- 
position) constraints (see Definition 12.11 below) . In this paper, we propose to use this reformulation 
leading to a new local optimization method for solving some classes of optimization problems involv- 
ing BMI constraints. We provide a practical algorithm and prove the convergence of the algorithm 
under certain standard assumptions. 

The algorithm proposed in this paper is very simple to implement by using available SDP software 
tools. Moreover, it does not require any globalization strategy such as line-search procedures to 
guarantee global convergence to a local minimum. The method still works in practice for nonsmooth 
optimization problems, where the objective function and the concave parts are only subdifferentiable, 
but not necessarily differentiable. Note that our method is different from the standard DCA approach 
in |28] since we work directly with positive semidefinite matrix inequality constraints instead of 
transforming into DC representations as in [2]. 

We show that our method is applicable to many control problems in static state/output feedback 
controller design. The numerical results are benchmarked using the data from COMPl c ib library. 
Note, however, that this method is also applicable to other nonconvex optimization problems with 
matrix inequality constraints which can be written as a convex-concave decomposition. 

Outline of the paper. The remainder of the paper is organized as follows. Section [5] provides 
some preliminary results which will be used in what follows. Section [3] presents the formulation of 
optimization problems involving convex-concave matrix inequality constraints and a fundamental 
assumption, Assumption AfT] The algorithm and its convergence results are presented in Section 0] 
Applications to control problems on static feedback controller design and numerical benchmarking 
are given in Section [5] The last section contains some concluding remarks. 

2 Preliminaries 

Let S p be the set of symmetric matrices of size p x p, <S?, and resp., <S? . be the set of symmetric 
positive semidefinite, resp., symmetric positive definite matrices. For given matrices X and Y in 



x = Ax + B\w + Bu, 
z = C\x + D n w + D 12 u, 
y = Cx + D 2 iw, 



(1.1) 




(1.2) 
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S p , the relation X h Y (resp., X <Y) means that I - 7 e 5j (resp., Y ~ X S+) and X >- Y 
(resp., 1 ^ 7) is 1 - 7 e (resp., F - X G <S£ + ). The quantity Ioy ; = trace(X T F) is an 
inner product of two matrices X and Y defined on S p , where trace(Z) is the trace of matrix Z. 

Definition 2.1. [33] A matrix- valued mapping G : R" — > S p is said to be positive semidcfinitc 
convex (psd-convex) on a convex subset G C 1" if for all t G [0, 1] and x, y G C, one has 

G(ta + (1 - t)y) r< tG(x) + (1 - i)G(y). (2.1) 

If (12.11) holds true for ^ instead of -< and i G (0, 1) then G is said to be strictly psd-convex on G. 
Alternatively, if we replace ^ in (|2.ip by >r then G is said to be psd-concave on G. It is obvious 
that any convex function / : M. n — >■ E is psd-convex (p = 1). 

A function / : R ra — > K is said to be strongly convex with the parameter p > if /(■) — ^pH ■ || 2 
is convex. 

The derivative of a matrix-valued mapping G at x is a linear mapping DG from K™ to M. pxp 
which is defined by 

" BC 

DG(x)h := V hi — (x), \/h G K™. 

For a given convex set X G K", the matrix- valued mapping G is said to be differentiable on 
a subset X if its derivative DG(x) exists at every x G X. The definitions of the second order 
derivatives of matrix- valued mappings can be found, e.g., in [25]. Let A : M. n —> S p be a linear 
mapping defined as Ax = 5Zi=i x iA-i, where Ai G 5 P for i = 1, . . . , n. The adjoint operator of A, 
A*, is defined as A* Z = {A x o Z, A 2 o Z, . . . , A n o Z) T for any Z G 5 P . 

Lemma 2.1. a) A matrix-valued mapping G is psd-convex on X if and only if for any v G W 
the function ip(x) := v T G(x)v is convex on X . 

b) A mapping G is psd-convex on X if and only if for all x and y in X, one has 

G(y)-G{x)yDG(x)(y-x). (2.2) 

Proof. The proof of the statement a) can be found in [29 . We prove b). Let <p(x) = v T G(x)v for 
any v G K p . If G is psd-convex then <p is convex. We have f{y) — f{x) > \7<p(x) T (y — x). Now, 
Vcp(x) T (y -x) = ££=i (ift - Xi)v T §£(x)v = v T [DG(x)(y - x)]v. Hence, v T \G{y) - G(x) - DG(x) 
(y — x)] v > for all v. We conclude that (|2.2[) holds. Conversely, if (|2.2[) holds then, for any v, we 
have v T \G{y) — G(x) — DG(x)(y — x)]v > 0, which is equivalent to ip(y) — tp(x) > \7tp(x) T (y — x). 
Thus ip is convex. By virtue of a) , the mapping G is psd-convex. □ 

For simplicity of discussion, throughout this paper, we assume that all the functions and matrix- 
valued mappings are twice differentiable on their domain |291 134 ] . However, this assumption can 
be reduced to the subdifferentiability of the objective function and the concave parts of the matrix- 
valued mappings. 

Definition 2.2. A matrix-valued mapping F : 1" — > S p is said to be a psd-convex-concave mapping 
if F can be represented as a difference of two psd-convex mappings, i.e. F(x) = G(x) — H(x), where 
G and H are psd-convex. The pair (H, G) is called a psd-DC (or psd-convex-concave) decomposition 
ofF. 

Note that each given psd-convex-concave mapping possesses many psd-convex-concave decom- 
positions. 

3 Optimization of convex-concave matrix inequality con- 
straints 

3.1 Psd-convex-concave decomposition of BMIs 

Instead of using the vector a; as a decision variable, we use from now on the matrix X as a matrix 
variable in K mx ". Note that any matrix X can be considered as an m x n-column vector by 
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vectorizing with respect to its columns, i.e. x = vec(X) := (Xu, X21, ■ ■ ■ , X mn ) T . The inverse 
mapping of vec is called mat. Since vec and mat are linear operators, the psd-convexity is still 
preserved under these operators. 

A mapping F : R pxq x S p -> S p given by F(X, Y) := XQ- l X T - Y, where Q G S q is symmetric 
positive definite, is called a Schur psd-conve^ mapping. 

Consider a bilinear matrix form 



F(X,Y) := X T Y + Y T X. 



(3.1) 



By using the Kronecker product, we can write F as vec(F(X,Y)) = (J„ <E> X T )vec(Y) + (I y <E> 
• XiUj), where I n , I y are appropriate identity matrices, ® denotes the Kronecker 
product. Hence, the vectorization of F(X, Y) is indeed a bilinear form of two vectors x := vec(X) 
and y :— vec(Y). 

The following lemma shows that the bilinear matrix form p. II) can be decomposed as a difference 
of two psd-convex mappings. 



Lemma 3.1. a) The mapping f(X) :— X T X . g(X) := XX T ) are psd-convex on 
mapping f(X) := X^ 1 is psd-convex on S+ + . 



The 



b) The bilinear matrix form X T Y + Y T X can be represented as a psd-convex- concave mapping 
of at least three forms: 



X T Y + Y T X = (X + Y) T {X + Y) - (X T X + Y T Y) 
= X T X + Y T Y - (X - Y) T (X - Y) 



(3.2) 
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[(X + Y) T (X + Y) - (X - Yf(X - Y)]. 



The statement b) provides at least three different explicit psd-convex-concave decompositions of 
the bilinear form X T Y + Y T X. Intuitively, we can see that the first decomposition has a "strong 
curvature" on the second term, while the second and the third decompositions have "less curvature" 
on the second term due in case of a compensation between X and Y. 

The following result will be used to transform Shur psd-convex constraints to LMI constraints. 



Lemma 3.2. 

to 



a) Suppose that A G S n . Then the matrix inequality BB T — A -< (^) is equivalent 



A B 
B T I 



b) Suppose that A 6 S n , D y 0, then we have: 

1 >■ (h) 



A - BB T C 
C T D 



y (b) 0. 



A 
B T 
C T 



(3.3) 



B C 
I O 
O D 



> (>=) 0. 



(3.4) 



The proof of this lemma immediately follows by applying Schur's complement and Lemma 12.11 
[S]. We omit the proof here. 

3.2 Optimization involving convex-concave matrix inequality constraints 

Let us consider the following optimization problem: 



min fix) 

s.t. Gi(x) - Hi(x) < 0, i 
x G fi, 



.,1, 



(3.5) 



1 Due to Schur's complement form 
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where / : R n — > R is convex, f2 C R™ is a nonempty, closed convex set, and Gj and Hi (i = 1, . . . , I) 
are psd-convex. Problem (|3.5p is referred to as a convex optimization with psd-convex-concave 
matrix inequality constraints. 

Let ft be a polyhedral in R™. Then, if / is nonlinear or one of the mappings Gj or Hi (i = 1, . . . , I) 
is nonlinear then (|3.5[) is a nonlinear semidefinite program. If Hi (i = 1, . . . , I) are linear then (|3.5[) 
is a convex nonlinear SDP problem. Otherwise, it is a nonconvex nonlinear SDP problem. 

Let us define L(x, A) := f(x) + Yli=i A» ° [Gi(x) — £/i(x)] as the Lagrange function of (13.51) . where 
Aj G <S P (i = 1, . . . , I) considered as Lagrange multipliers. The generalized KKT condition of (|3.5[) 
is presented as: 

e V/(x) + ELi £>[G<(x) - H^x)]*^ + N a (x), 
Gi(x)-Hi(x) 1 0, Ai h 0, (3.6) 
[G i (x)-H i (x)]oA i =Q, i = l,...,l. 



Here, ATn(x) is the normal cone of fi at x defined as 

N n (x) := 



_J {u; e R™ | w T (y-x) > 0, Vy G 0} , ifxGO, 
0, otherwise. 

A pair (x*,A*) satisfying (|3.6p is called a KKT point, x* is called a stationary point and A* is the 
corresponding multiplier of Q3.5p . The generalized optimality condition for nonlinear semidefinite 
programming can be found in the literature (e.g., \TZ\ ). 
Let us denote by 

V:={xen\ G l {x) - Hi(x) < 0, 1 = 1,..., I}, (3.7) 
the feasible set of f|3.5[) and ri(D) is the relative interior of V which is defined by 

ri(D) := {x G ri(fi) | G l {x) - H,(x) -< 0, i = 1, . . . , m} , 

where ri(fi) is the set of classical relative interiors of fi [51 [TH]. The following condition is a funda- 
mental assumption in this paper. 

Assumption A.l. ri(D) is nonempty. 

Note that this assumption is crucial for our method, because, as we shall see, it requires a strictly 
feasible starting point x° G ri(2?). Finding such a point is in principle not an easy task. However, 
in many problems, this assumption is always satisfied. In Section [S] we will propose techniques to 
determine a starting point for the control problems under consideration. 



4 The algorithm and its convergence 

In this section, a local optimization method for finding a stationary point of problem (13.51) is pro- 
posed. Motivated from the DC programming algorithm developed in [2 8) and the convex-concave 
procedure in |30j for scalar functions, we develop an iterative procedure for finding a stationary 
point of p. 51) . The main idea is to linearize the nonconvex part of the psd-convex-concave matrix 
inequality constraints and then transform the linearized subproblem into a quadratic semidefinite 
programming problem. The subproblem can be either directly solved by means of interior point 
methods or transformed into a quadratic problem with LMI constraints. In the latter case, the 
resulting problem can be solved by available software tools such as Sedumi [21] and SDPT3 [3"3"] . 

4.1 The algorithm 

Suppose that x fe G ft is a given point, the linearized problem of (|3.5p around x k is written as 
f mm {f k ( x ):=f( x ) + f\\Q k ( x - Xk )f 2 } 

< s*t. G l {x)-H l (x k )~DH l (x k )(x-x k )d : 0,i = l,...,l, (4.1) 
I x G f2. 
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Here, we add a regularization term into the objective function of the original problem, where Qk 
is a given matrix that projects x — x^ in a certain subspace of MP and pk > is a regularization 
parameter. Since Gi (i = 1, . . . , /) are psd-convex and the objective function is convex, problem (|4.1j) 
is convex. The linearized convex- concave SDP algorithm for solving (|3.5p is described as follows. 

Algorithm 1. 

Initialization: Choose a positive number po and a matrix Qq G 5™. Find an initial point x° G ri(2?). 
Set fc := 0. 

Iteration fc: For fc = 0, 1, . . . . Perform the following steps: 

Step 1: Solve the convex semidefmite program (|4.1[) to obtain a solution x k+1 . 

Step 2: If ||a; +1 — x k \\ < e for a given tolerance e > then terminate. Otherwise, update pk 
and Qk (if necessary), set fc := fc + 1 and go back to Step 1. 

The following main property of the method makes an implementation very easy. If the initial 
point x a belongs to the relative interior of the feasible set T>, i.e. x° G ri(2?), then Algorithm [1] 
generates a sequence x k which still belongs to T>. Consequently, no line-search procedure is needed 
to ensure the global convergence. 

This property follows from the fact that the linearization of the concave part —Hi is its an upper 
approximation of this mapping (in the sense of positive semidefmite cone), i.e. 

-Hi(x) < -Hi{x k ) - DHi(x k )(x - x k ),Vx G O, 

which is equivalent to 

Gi(x) - Hi(x) r< Gi(x) - Hi{x k ) - DHi(x k )(x - x k ),Vx G fi. 

Hence, if the subproblem (I4.ip has a solution x k+1 then it is feasible to f|3 . 5f) . Geometrically, 
Algorithm [1] can be seen as an inner approximate method. 

The main tasks of an implementation of Algorithm [1] consist of: 

• determining an initial point a; G ri(X>), and 

• solving the convex semidefmite program (I4.1[) repeatedly. 

As mentioned before, since T> is nonconvex, finding an initial point x° in ri(Z?) is, in principle, not an 
easy task. However, in some practical problems, this can be done by exploiting the special structure 
of the problem (see the examples in Section [5J . 

To solve the convex subproblem (|4.1[) , we can either implement an interior point method and 
exploit the structure of the problem or transform it into a standard SDP problem and then make 
use of available software tools for SDP. The regularization parameter pk and the projection matrix 
Qk can be fixed at appropriate choices for all iterations, or adaptively updated. 

Lemma 4.1. If x k is a solution of (|4.ip linearized at x k then it is a stationary point of (|3.5[) . 

Proof. Suppose that A fe+1 is a multiplier associated with x k , substituting x k into the generalized 
KKT condition (TTJ) of (|4~Tj) we obtain (gig) . Thus x k is a stationary point of (|3~5|) . □ 

4.2 Convergence analysis 

In this subsection, we restrict our discussion to the following special case. 

Assumption A. 2. The mappings £?, (i = 1, . . . , I) are Schur psd-convex and f2 is formed by a finite 
number of LMIs. In addition, / is convex quadratic on R n with a convexity parameter pf > 0. 

This assumption is only technical for our implementation. If the mapping Gi is Schur psd-convex 
then the linearized constraints of problem (|4.1I) can directly be transformed into LMI constraints 
(see Lemma l3.2j) . In practice, G; (i = 1, . . . , I) can be a general psd-convex mappings and / can be 
a general convex function. 
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Under Assumption the convex subproblem (|4.1I) can be transformed equivalently into a 
quadratic semidefinite program of the form: 




\z T Bz + h T z 

(4.2) 

A(z) + C^0, 

where A is a linear mapping from R™ 2 to S p * , C G S Pz and B is a symmetric matrix, by means of 
Lemma 13.21 

A vector z is said to satisfy the Slater condition of (|4.2j) if A(z) + C -< 0. Suppose that the triple 
(z, V, S) satisfies the KKT condition of (|4. 2[) (see [H]), where z is a primal stationary point, V is a 
Lagrange multiplier and 5 is a slack variable associated with z and V. Then, problem (|4.2[) is said 
to satisfy the strict complementarity condition at (z, V , S) if V + S y 0. 

Let 2 be a stationary point of (]4.2j) . We say that ^ p 6 R™ z is a feasible direction to (|4.2p if 
2 + ep is a feasible point of (|4.2I) for all e > sufficiently small. As in [TT| , we assume that the second 
order sufficient condition holds for (|4.2[) at z with modulus /i > if for all feasible directions p at 
z with p T (h + Bz) = 0, one has p 7 Bp > /i||p|| 2 . We say that the convex problem (|4.2p is solvable 
and satisfies the strong second order sufficient condition if there exists a KKT point (z, V, S) of the 
KKT system of (|4.2[) satisfies the second order sufficient condition and the strict complementary 
condition. 

Assumption A. 3. The convex subproblem (|4.f p is solvable and satisfies the strong second order 
sufficient condition. 

Assumption is standard in optimization and is usually used to investigate the convergence of 
the algorithms [TD1 QH H3 . 

The following lemma shows that Ax k :— x k+1 — x k is a descent direction of problem (|3.5[) whose 
proof is given in the Appendix. [7] 

Lemma 4.2. Suppose that {(x k , A k )}k>o is a sequence generated by Algorithm^ Then: 

a) The following inequality holds for k > 0: 

f{x k+l ) - f(x k ) < -^\\x k+1 x k \\ 2 2 Pk \\Q k (x k+1 - x k )\\l (4.3) 
where Pf is the convexity parameter of f. 

b) If there exists at least one constraint i$, io S {1,2,...,/}, to be strictly feasible at x k , i.e. 
G l0 (x k ) - H lQ (x k ) -< 0, then f{x k+1 ) < f(x k ) provided that A k+1 y 0. 



c) If pk > and Qk is full-row-rank then Ax k is a sufficiently descent direction of (|3.5I) . 

The following theorem shows the convergence of Algorithm [T] in a particular case. 

Theorem 4.1. Under Assumptions ^Ql and suppose that f is bounded from below on T>, 
where T> is assumed to be bounded in R". Let {(a^A^)} be a sequence generated by Algorithm 
[7] starting from x° G ri(T>). Then if either f is strongly convex or pk = p > and Qk = Q is 
full-row-rank for all k > then every accumulation point (x*,A*) of {(x k ,A k )} is a KKT point of 
p.5p . Moreover, if the set of the KKT points of p.5p is finite then the whole sequence {(x k ,A k )} 
converges to a KKT point of Q3.5P . 



Proof. Let M{x°) := {x k } be a sequence of the sample points generated by Algorithm [T] starting 
from x . For a given x G fl, let us define the following mapping: 

A sol (x) :=argmin{ f(y) + P - \\Q(y - x) \\ 2 2 \ y G SI, (4.4) 
Gi(y) - Hi(x) - DHi(x)(y - x) < 0, i = 1, . . . ,m}. 
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Then, A so \ is a multivalued mapping and it can be considered as the solution mapping of the convex 
subproblem (|4.1[) . Note that the sequence {x k } generated by Algorithm [T] satisfies x k+1 £ ^4 so i (x* ) . 
We first prove that A so i is a closed mapping. Indeed, since the convex subproblem (|4.1[) satisfies 
Slater's condition and has a solution that satisfies the strict complementarity and the second order 
sufficient condition, applying Theorem 1 in |Ilj we conclude that the mapping ^4 so i is differentiable 
in a neighborhood of the solution. Consequently, it is closed due to the compactness of T>. 

On the other hand, since / is either strongly convex or p^ = p > for all k > and Qk = Q 
is full-row-rank, it follows from Lemma 14.21 that the objective function / is strictly monotone on 
M(x°). Since M(x°) C D and V is compact, M(x°) is also compact. Applying Theorem 2 in [53] 
we conclude that every limit points of the sequence {x k } belongs to the set of stationary points S* . 
Moreover, S* is connected and if S* is finite then the whole sequence {x k } converges to x* in S*. □ 

Remark 4.1. The condition that / is quadratic in Assumption [2] can be relaxed to / being twice 
continuously differentiable. However, in this case, we need a direct proof for Theorem l4. II instead of 
applying Theorem 1 in 

5 Applications to robust controller design 

In this section, we apply the method developed in the previous section to the following static 
state/output feedback controller design problems: 

1. Sparse linear static output feedback controller design; 

2. Spectral abscissa and pseudo-spectral abscissa optimization; 

3. 'Hi optimization; 

4. "Hoo optimization; 

5. and mixed 7^2 /%oo synthesis. 

We used the system data from [21 [37] and the COMPLib library [2Q]. All the implemen- 
tations are done in Matlab 7.11.0 (R2010b) running on a PC Desktop Intel(R) Core(TM)2 
Quad CPU Q6600 with 2.4GHz and 3Gb RAM. We use the YALMIP package [22] with the 
SeDuMi 1.1 solver [3T] to solve the LMI optimization problems arising in Algorithm [1] at 
the initial phase (Phase 1) and subproblem (|4.1[) . The Matlab codes can be downloaded at 
http://www.kuleuven.be/optec/software/BMIsolver, We also benchmarked our method with 
various examples and compared our results with HIFOO [T3] and PENBMI [TS] for all control prob- 
lems. HIFOO is an open-source Matlab package for fixed-order controller design. It computes a 
fixed-order controller using a hybrid algorithm for nonsmooth, nonconvex optimization based on 
quasi-Newton updating and gradient sampling. PENBMI 15 is a commercial software for solving 
optimization problems with quadratic objective and BMI constraints. PENBMI is free licensed for 
academic purposes. We initialized the initial controller for HIFOO and the BMI parameters for 
PENBMI to the initial values of our method. As we shall see, we can reformulate the spectral 
abscissa optimization problem as a rank constrained LMI problem. Therefore, we also compared 
our results with LMIRank [Mj , a MATLAB toolbox for solving rank constrained LMI problems, for 
the spectral abscissa optimization. 

Note that all problems addressed here lead to at least one BMI constraint. To apply the method 
developed in the previous section, we propose a unified scheme to treat these problems. 

Scheme A.l. 

Step 1. Find a convex-concave decomposition of the BMI constraints as G(x) — H(x) -< 0. 
Step 2. Find a starting point x° e ri(X>). 

Step 3. For a given x , linearize the concave part to obtain the convex constraint G(x) — 
Hk{x) di 0, where is the linearization of H at x k . 
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Step 4- Reformulate the convex constraint as LMI constraint by means of Lemma l3.2l 
Step 5. Apply Algorithm [1] with SDP solver to solve the problem. 



5.1 Sparse linear constant output-feedback design 

Let us consider a BMI optimization problem of sparse linear constant output-feedback design given 
as: 

min -on + Ym-i Wv I 

s.t (A + BFC) T P + P(A + BFC) + 2aP -< 0, P = P T , P y 0. 

Here, matrices A, B, C are given with appropriate dimensions, P and F are referred as variables 
and a > is a weighting parameter. The objective function consists of two terms: the first term 
aa is to stabilize the system (or to maximize the decay rate) and the second one is to ensure the 
sparsity of the gain matrix F. This problem is a modification of the first example in [14]. Let us 
illustrate Scheme A{T]for solving this problem. 

Step 1: Let Bp := A + BFC + al, where I is the identity matrix. Then, applying Lemma 13. II 
we can write 

(A+BFC) T P+P(A+BFC) + 2aP = B T F P + PB F 

= bIb f + P t P-{B f -P) t (B f -P), (5.2) 
- \ [{Bf + P) T (B F + P) - {B F - P) T (B F - P)] . (5.3) 



In our implementation, we use the decomposition f|5 .3[) . 
If we denote by 

G(a,P,F) :=^(B F + P) T (B F + P), and H(a, P, F) := ^(B F - P) T (B F - P), (5.4) 



then the BMI constraint in (15.11) can be written equivalently as a psd-convex-concave matrix inequal- 
ity constraint (of a variable x formed from (a,P,F) as x := (a, vec(P) T , vec(F) T ) T ) as follows: 

G{a,P,F) - H{a,P,F) ^ 0. (5.5) 

Note that the objective function of (|5.ip is convex but nonsmooth which is not directly suitable 
for the SSDP approach in [9], but, the nonconvex problem (|5.ip can be reformulated in the form of 
p.5[) using slack variables. 

Steps 2-5: The implementation is carried out as follows: 

Phase 1. {Determine a starting point a; G ri(T>)). Set F° := 0, a := — cto (A T + A)/2 where 
ao(') is the maximum real part of the eigenvalues of the matrix, and compute P = P° as the 
solution of the LMI feasibility problem 

(A + BF°C) T P + P(A + BF°C) + 2a°P -< 0. (5.6) 

The above choice for (a°,F°) originates from the property that P° = I renders the left hand 
size of (|5.6[) negative semi-definite (but not negative definite). 

Phase 2. Perform Algorithm [T] with a starting point x° found at Phase 1. 

Let us now illustrate Step 4 of Scheme AfT] After linearizing the concave part of the convex-concave 
reformulation of the last BMI constraint in (I5.1[) at (F k , P k ,a k ) we obtain the linearization: 

(A+BFC+aI+P) T {A+BFC+aI+P) - H k {F. P, a) -< 0, (5.7) 

where Hk(F,P,a) is a linear mapping of F, P and a. Now, applying Lemma [3.21 (|5.7p can be 
transformed to an LMI constraint: 



H k (F,P,a) {A+BFC+aI+P) T ' 
(A + BFC+al+P) I 



y 0. 
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With the above approach we solved problem (|5.ip for the same system data as in [13]. Here, matrices 
A, B and C are given as: 



A= 



2.45 -0.90 
-0.12 -0.44 

2.07 -1.20 
-0.59 0.07 
-0.74 -0.23 



1.53 
-0.01 
-1.14 

2.91 
-1.19 



-1.26 
0.69 
2.04 
-4.63 
-0.06 



1.76 
0.90 
-0.76 
-1.15 
-2.52 



B= 



0.81 
-0.34 
-1.32 
-2.11 

0.31 



-0.79 
-0.50 
1.55 
0.32 
-0.19 



0.00 
0.06 

-1.22 
0.00 

-1.09 



0.00 
0.22 
-0.77 
-0.83 
0.00 



-0.95 
0.92 

-1.14 
0.59 
0.00 



C= 



0.00 0.00 0.16 0.00 -1.781 

1.23 -0.38 0.75 -0.38 -0.00 

0.46 0.00 -0.05 0.00 0.00 

0.00 -0.12 0.23 -0.12 1.14 



The weighting parameter a is chosen by a = 3. Algorithm [T] is terminated if one of the following 
conditions is satisfied: 

• subproblcm (14. ip encounters a numerical problem; 



|Aa 



I) < 10 



-3. 



• the maximum number of iterations, JsT ma x, reaches; 

• or the objective function is not significantly improved after two successive iterations (i.e. 
|jfc+i _ jfc| < io-4(i + lyfe^ f or some k = k&ndk = k + l, where f k := f(x k )). 

In this example, Algorithm [T] is terminated after 15 iterations, whereas the objective function is not 
significantly improved. However, after the 2 nd iteration, matrix F only has 3 nonzero elements, while 
the decay rate a is 1.17316. This value is much higher than the one reported in [T3], a — 0.3543 
after 6 iterations. We obtain the gain matrix F as 



F = 



0.6540 0.0000 0.0000 0.0000 

0.0000 -0.4872 0.0000 0.0000 

0.0000 0.0000 0.0000 0.0000 

0.0000 0.0000 0.0000 0.0000 

0.0000 1.1280 0.0000 0.0000 



With this matrix, the maximum real part of the eigenvalues of the closed- loop matrix in (11.21) . 
Ap := A + BFC, is a {A F ) := -1.40706. Simultaneously, a {A T F P + PA F + 2aP) = -0.327258 < 
and ao(P) = 0.587574 > 0. Note that ao{Ap) ^ —a due to the in-activeness of the BMI constraint 
in (|5.1I) at the 2nd iteration step. 



5.2 Spectral abscissa and pseudo-spectral abscissa optimization 

One popular problem in control theory is to optimize the spectral abscissa of the closed-loop system 
x = (A + BFC)x. Briefly, this problem is presented as an unconstrained optimization problem of 
the form: 

min a Q (A + BFC), (5.8) 

where a Q (A + BFC) := sup {Re(A) | A E A(A + BFC)} is the spectral abscissa of A + BFC, Re(A) 
denotes the real part of A e C and A(A + BFC) is the spectrum of A + BFC. Problem (pT8)) 
has many drawbacks in terms of numerical solution due to the nonsmoothness and non-Lipschitz 
continuity of the objective function ao [7]. 

In order to apply the method developed in this paper, we reformulate problem (|5.8[) as an 
optimization problem with BMI constraints [JJ [3T] : 

{max B 
P > F >P (5.9) 
s.t. (A + BFC) P + P(A + BFC) + 2BP -< 0, P = P T , P >- 0. 

Here, matrices A e W nxn , B 6 M" xn ", C 6 M"« x " are given. Matrices P £ R nx ™ and F G M"" x ™« 
and the scalar j3 are considered as variables. If the optimal value of (|5.9[) is strictly positive then 
the closed-loop feedback controller u = Fy stabilizes the linear system x = {A + BFC)x. 

Problem (|5.9p is very similar to (|5.ip . Therefore, using the same trick as in (15.11) . we can 
reformulate (|5.9p in the form of (I3.5p . More precisely, if we define Bp := A + BFC + (31 then the 
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bilinear matrix mapping Ap^P + PAp can be represented as a psd-convex-concave decomposition of 
the form (|5.3p and problem (I5.9[) can be rewritten in the form of (|3.5I) . We implement Algorithm 
[T] for solving this resulting problem using the same parameters and the stopping criterions as in 
Subsection 15. II In addition, we regularize the objective function by adding the term ^-\\F — F k \\p + 
^-\\P — P k \\p, with pf = pp = 1CP 2 . The maximum number of iterations -ftT ma x is set to 150. 

We test for several problems in COMPl e ib and compare our results with the ones reported by 
HIFOO, PENBMI and LMIRank. For LMIRank, we implement the algorithm proposed in [25] . 
We initialize the value of the decay rate a at 10~ 4 and perform an iterative loop to increase a as 
a k+1 := a + 0.1. The algorithm is terminated if either the problems (12) or (21) in [55] with a 
correspondence a can not be solved or the maximum number of iterations K max = 100 is reached. 
The numerical results of four algorithms are reported in Table [TJ Here, we initialize the algorithm 



Table 1: Computational results for (j53)l in COMPLib 



Problem 


Other Results, ao(Ap) 


Results and Performances 


Name 


otoiA) 


HIFOO 


LMIRANK 


PENBMI 


a (A F ) 


Iter 


time [s] 


AC I 


0.00(1 


-0.2061 


-8.4766 


-7.0758 


-0.8535 


41 


12.44 


81 A(J4 


2.579 


-0.0500 


-0.0500 


-0.0500 


-0.0500 


14 


4.60 


AC5" 


0.999 


-0.7746 


-1.8001 


-2.0438 


-0.7389 


28 


63.33 


AC 7 


0.172 


-0.0322 


-0.0204 


0.0896 


-0.0673 


150 


111.46 


'2 ACS 


0.012 


-0.1968 


-0.4447 


0.4447 


-0.0755 


24 


21.95 


6 acs) 


0.012 


-0.3389 


-0.5230 


-0.4450 


-0.3256 


78 


74.57 


7AC11 


5.451 


-0.0003 


-5.0577 




-3.0244 


61 


38.44 


52AC12 


0.580 


-10.8645 


-9.9658 


-1.8757 


-0.3414 


150 


86.72 


100 HE1 


o/m 


-0.2457 


-0.2071 


-0.2468 


-0.2202 


150 


87.64 


32.8000. 4HE3 


0.087 


-0.4621 


-2.3009 


-0.4063 


-0.8702 


47 


48.80 


700.8300, 25HE4 


0.234 


-0.7446 


-1.9221 


-0.0909 


-0.8647 


63 


71.66 


21HE5 


0.234 


-0.1823 




-0.2932 


-0.0587 


150 


178.71 


HE6 


0.234 


-0.0050 


-0.0050 


-0.0050 


-0.0050 


12 


41.00 


2 REM 


1.991 


-16.3918 


-5.9736 


-1.7984 


-3.8599 


77 


79.23 


61 REA2 


2.011 


-7.0152 


-10.0292 


-3.5928 


-2.1778 


40 


39.15 


REA3 


0.000 


-0.0207 


-0.0207 


-0.0207 


-0.0207 


150 


362.21 


D1S2 


1.675 


-6.8510 


-10.1207 


-8.3289 


-8.4540 


28 


37.28 


D1S4 


1.442 


-36.7203 


-0.5420 


-92.2842 


-8.0989 


72 


124.23 


WEC1 


0.008 


-8.9927 


-8.7350 


-0.9657 


-0.8779 


150 


305.36 


1H 


0.000 


-0.5000 


-0.5000 


-0.5000 


-0.5000 


7 


39.41 


CSE1 


0.000 


-0.4509 


-0.4844 


-0.4490 


-0.2360 


38 


158.67 


TF1 


0.000 






-0.0618 


-0.1544 


56 


137.98 


TF2 


0.000 






-l.Oe-5 


-l.Oe-5 


8 


20.41 


TF3 


0.000 






-0.0032 


-0.0031 


93 


237.93 


NN1 


3.600' 


-3.0458 


-4.4021 


-4.3358 


-0.8746 


12 


37.53 


53.1700. 65JNJN5" 


0.420 


-0.0942 


-0.0057 


-0.0942 


-0.0913 


11 


42.62 


NN9 


3.281 


-2.0789 


-0.7048 




-0.0279 


33 


111.41 


NN13 


1.945 


-a. 251a 


-4.5810 


-9.0741 


-3.4318 


150 


572.50 


NN15 


0.000 


-6.9983 


-11.0743 


-0.0278 


-0.8353 


150 


524.80 


NJN17 


1.170 


-0.6110 


-U.5130 




-0.6008 


99 


342.67 



in HIFOO with the same initial guess F° — 0. Since PENBMI and our methods solve the same BMI 
problems, they are initialized by the same initial values for P, F and /3. 

The notation in Table [T] consists of: Name is the name of problems, ao(A), ao(Ap) are the 
maximum real part of the eigenvalues of the open-loop and closed-loop matrices A, Ap, respectively; 
iter is the number of iterations, time [s] is the CPU time in second. The columns titled HIFOO, 
LMIRank and PENBMI give the maximum real part of the eigenvalues of the closed-loop system for 
a static output feedback controller computed by available software HIFOO [13], LMIRank [55] and 
PENBMI [15], respectively. Our results can be found in the sixth column. The entries with a dash 
sign indicate that there is no feasible solution found. Algorithm 1 fails or makes only slow progress 
toward to a local solution with 6 problems: AC18, DIS5, PAS, NN6, NN7, NN12 in COMPl e ib. 
Problems AC 5 and NN5 are initialized with a different matrix F° to avoid numerical problems. 

Note that Algorithm[T]as well as the algorithms implemented in HIFOO, LMIRank and PENBMI 
are local optimization methods, which only report a local minimizcr and these solutions may not be 
the same. To apply the LMIRank package for solving problem (15.91) . we have used a direct search 
procedure for finding a. The computational time of this procedure is very high compared with the 
other methods. 

To conclude this subsection, we show that our method can be applied to solve the optimization 
problem of pseudo-spectral abscissa in static feedback controller designs. This problem is described 
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as follows (see [7HH]): 



max 8 

~2/3P+A' p P+PA F +fiI-ujI z eP 
eP ujI 



s.t 



XO, P >- 0, P = P T , /i > 0, 



(5.10) 



where Ap = A + BFC as before and w < 0. 

Using the same notation Pp = A + BFC + BI as in (|5.9[) and applying the statement b) of 
Lemma 13.21 the BMI constraint in this problem can be transformed into a psd-convex-concave one: 



BpB F + P J P+( f i-uj)I eP 
eP w7 



(B F -P)\B F -P) 




-<0. 



If we denote the linearization of (Bp — P) T (B F — P) at the iteration fc by H k , i.e. Pfc = (Pp — 
P) T (B% - P fc ) + (B F - P k ) T (B F - P) - (P| - P k ) T (B F - P fe ), then the linearized constraint in 
the subproblem (|4.ip can be represented as an LMI thanks to Lemma \'3. 2 1 



H k + (uj- 




p 


-eP 


B F 


I 








P 





/ 





-sP 








-LOI 



y o. 



Hence, Algorithm [T] can be applied to solve problem (15 . 10|) . 

Remark 5.1. If we define F := BFC then the bilinear matrix mapping A^P + PAp can be rewritten 



as 



A F P 



PAp = - [(P + F) T (P + P) - (P - F) T (P - F)] - A T P - PA. 



Using this decomposition, one can avoid the contribution of matrix A on the bilinear term. Conse- 
quently, Algorithm [1] may work better in some specific problems. 



5.3 TL2 optimization: BMI formulation 

In this subsection, we consider an optimization problem arising in the T-ii synthesis of the linear 
system (ll.lj) . Let us assume that Dyx = and P21 = 0, then this problem is formulated as the 
following optimization problem with BMI constraints [2D] , 



mm 

F.Q.X 



trace(A) 

s.t. (A + BFC)Q + Q(A + BFC) T + B X B\ < 0, 



X CiQ 
QCf Q 



h0, QyO. 



(5.11) 



Here, we also assume that BiBf is positive definite. Otherwise, we use BiBf + el instead of PiPf 
with e = 1(T 5 in (pUTi 

In order to apply Algorithm [T] for solving problem (|5.11l) . a starting point x° G ri(P) is required. 
This task can be done by performing some extra steps called Phase 1. The algorithm is now split 
in two phases as follows. 
Phase 1: (Determine a starting point x°). 

Step 1. If a (A + A T ) < then we set P° := 0. Otherwise, go to Step 3. 
Step 2. Solve the following optimization problem with LMI constraints: 



mintrace(A') s.t. A F oQ + QAp + B X B{ -< 0, 



X CiQ 
QCf Q 



>-0, Q^0, (5.12) 



where A F o := A + BF°C. If this problem has a solution Q° and X° then terminate Phase 1 
and using P° together with Q°,X° as a starting point x° for Phase 2. Otherwise, go to Step 
3. 
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Step 3. Solve the following feasibility problem with LMI constraints: 



Find P >~ and K such that: 



PA + A T P + KC + C T K T PB X 
BfP -^I, 



< 0, 



C[ p 



x d 



h0, 



to obtain K* and P* , where <7o is a given regularization factor. Compute F* := B + (P*)~ 1 K*, 
where B + is a pseudo-inverse of B, and resolve problem (|5.12j) with F° := F*. If problem 
(|5.12[) has a solution Q° and X° then terminate Phase 1 and set x° := (F°, Q°, X°). Otherwise, 
perform Step 4. 

Step 4- Apply the method in Subsection 15.21 to solve the following BMI feasibility problem: 

Find F and QyO such that : (A + BFC)Q + Q(A + BFC) T + BiBj -< 0. (5.13) 

If this problem has a solution F° then go back to Step 2. Otherwise, declare that no strictly 
feasible point is found. 

Phase 2: (Solve problem (|5.11[) V Perform Algorithm [JJ with the starting point a; found at Phase 



Note that Step 3 of Phase 1 corresponds to determining a full state feedback controller and 
approximating it subsequently with an output feedback controller. Step 4 of Phase 1 is usually 
expensive. Therefore, in our numerical implementation, we terminate Step 4 after finding a point 
such that a ((A + BFC)Q + Q(A + BFC) T + B x Bj) < -0.1. 

Remark 5.2. The algorithm described in Phase 1 is finite. It terminates either at Step 4 if no feasible 
point is found or at Step 2 if a feasible point is found. Indeed, if a feasible matrix F° is found at Step 
4, the first BMI constraint of (|5.12l) is feasible with some Q y 0. Thus we can find an appropriate 
matrix X such that X — CQC T >- 0, which implies the second LMI constraint of (I5.12j) is satisfied. 
Consequently, problem (|5.12[) has a solution. 

The method used in Phase 1 is heuristic. It can be improved when we apply to a certain problem. 
However, as we can see in the numerical results, it performs quite acceptable for majority of the test 
problems. 

In the following numerical examples, we implement Phase 1 and Phase 2 of the algorithm using 
the decomposition 



for the BMI form at the left-top corner of the first constraint in (|5.1ip . The regularization parameters 
and the stopping criterion for Algorithm [JJ are chosen as in Subsection 15 . 1 1 and K max — 300. We test 
the algorithm for many problems in COMPl c ib and the computational results are reported in Table 
[2l For the comparison purpose, we also carry out the test with HIFOO [13, and PENBMI [15], and 
the results are put in the columns marked by HIFOO and PENBMI in Table [51 respectively. The 
initial controller for HIFOO is set to F° and the BMI parameters for PENBMI are initialized with 
(F, Q, X) = (F°, Q°, X°). Here, n, n y , n z , n w , n u are the dimensions of problems, the columns titled 
HIFOO and PENBMI give the H2 norm of the closed-loop system for the static output feedback 
controller computed by HIFOO and PENBMI; iter and time [s] are the number of iterations and 
CPU time in second of Algorithm [TJ respectively, included Phase 1 and Phase 2. Problems marked 
by "b" mean that Step 4 in Phase 1 is performed. In Table [5J we only report the problems that were 
solved by Algorithm [TJ The numerical results allow us to conclude that Algorithm [TJ PENBMI and 
HIFOO report similar values for majority of the test problems in COMPl e ib. 

If D12 ^ then the second LMI constraint of (|5. lip becomes a BMI constraint: 



1. 



□ 




X 

Q(C 1 + D 12 FC) T 



(Ci + D 12 FC)Q 

Q 



bo, 



(5.14) 
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Tabic 2: H2 synthesis benchmarks on COMPl ib plants 



Problem 


Other Hcsults, H.2 


Results and Fcrrormanccs 




Tlx 










HlrOO 






lter 


t imeLs J 




5 


3 


3 


2 


3 


0.0250 


0.0061 




3 


3.380 


AC2 6 


5 


3 


3 


5 


3 


0.0257 


0.0075 


0.0540 


3 


3.390 


— A i 


— c — I 

o 


— ' 

4 


z 


— f — i 




— f — 








TTTTV 

Z.i. i. i. t 




7'i '-i«0 


AC4 


4 




I 


2 


2 


1 1 096Q 




1 1 096Q 
J. .1. . \J Z u u 


2 


990 


A(J6 


7 


4 


2 


7 


7 


2.8648 


2.8648 


2.8664 


153 


124.230 


AC 7 


9 


2 


1 


1 


4 


0.0172 


0.0162 


0.0176 


1 


0.670 


AC8 


9 


5 


1 


2 


10 


0.6330 


0.7403 


0.6395 


300 


282.760 


AC12" 


4 


4 


3 


1 


3 


0.0022 


0.0106 


0.0992 


36 


28.540 


AC15" 


4 


3 


2 


6 


4 


1.5458 


1.4811 


1.5181 


264 


85.390 


ACie" 


4 


4 


2 


6 


4 


1.4769 


1.4016 


1.4427 


300 


99.790 


AC17 


4 


2 


1 


4 


4 


1.5364 


1.5347 


1.5507 


171 


49.350 


HE2 


4 


2 


2 


4 


4 


3.4362 


3.4362 


4.7406 


263 


97.310 


HE3" 


8 


6 


4 


10 


1 


0.0197 


0.0071 


0.1596 


249 


217.360 


HE4" 


8 


6 


4 


12 


8 


6.6436 


6.5785 


7.1242 


300 


412.830 


HEA1 


4 


3 


2 


4 


4 


0.9442 


0.9422 


1.0622 


249 


80.810 


REA2" 


4 


2 


2 


4 


4 


1.0339 


1.0229 


1.1989 


300 


101.730 


D1S1 


8 


4 


4 


8 


1 


0.6705 


0.1174 


0.7427 


300 


255.810 


U1S2 


3 


2 


2 


3 


3 


0.4013 


0.3700 


0.3819 


4 


1.370 


D1S3 


6 


4 


4 


6 


6 


0.9527 


0.9434 


1.0322 


300 


210.470 


U1S4 


e 


6 


4 


6 


6 


1.0117 


0.9696 


1.0276 


300 


210.690 


WECl" 


10 


4 


3 


10 


10 


7.3940 


8.1032 


12.9093 


119 


190.150 


WEC2 6 


10 


4 


3 


10 


10 


6.7908 


7.6502 


12.2102 


261 


407.470 


AGS 


12 


2 


2 


12 


12 


6.9737 


6.9737 


6.9838 


18 


28.830 


BUT 1 


11 


3 


3 


6 


1 


0.0024 




0.0017 


51 


64.410 


MFP 


4 




3 


4 


4 


6.9724 


6.9724 


7.0354 


300 


114.560 


PSM 


7 


3 


2 


5 


2 


0.0330 


0.0007 


0.1753 


300 


217.250 


EB2 6 


10 


1 


1 


2 


2 


0.0640 


0.0084 


0.1604 


114 


131.380 


EB3 


10 


1 


1 




2 


0.0732 


0.0072 


0.0079 


7 


6.240 


TFl" 


7 


4 


2 


4 


1 


0.0945 




0.1599 


192 


166.810 


TF2 


7 


3 


2 


4 


1 


11.1803 




11.1803 


3 


2.810 


TF3" 




3 


2 


4 


1 


0.1943 


0.1424 


0.2565 


138 


128.250 


NN2 


2 


1 


i 


2 


2 


1.1892 


1.1892 


1.1892 


4 


1.090 


NN4 


4 


3 


2 


4 


4 


1.8341 


1.8335 


1.8590 


222 


67.260 


NN8 


3 


2 


2 


3 


3 


1.5152 


1.5117 


1.5725 


183 


50.170 


NN11 


16 


5 


3 


3 


3 


0.1178 


0.0790 


0.1263 


39 


91.070 


NN13 6 


6 


2 


2 


3 


3 


26.1012 


26.1314 


62.3995 


138 


112.750 


NN14" 


6 


2 


2 


3 


3 


26.1448 


26.1314 


62.3995 


138 


110.020 


NN15 


3 


2 


2 


4 


1 


0.0245 




0.0210 


6 


2.060 


NN16" 


8 


4 


4 


4 


8 


0.1195 


0.1195 


0.1195 


3 


23.030 


NN17 


3 


1 


2 


2 


1 


3.2530 


3.2404 


3.3329 


300 


88.770 
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which is equivalent to X — CfQC f >: 0, where Cf '■— C\ + D\ 2 FC. Since f(Q) :— Q^ 1 is convex 
on S" ; !|_ (see Lemma 13. II a)), this BMI constraint can be reformulated as a convex-concave matrix 
inequality constraint of the form: 



' X 


C F ' 




O ' 


r T 
y F 





+ 


Q-\ 



By linearizing the concave term -f(Q) at Q = Q k as (QT 1 - (Q*) _1 (Q - Q fe )( ( 9 fc )" 1 (see [6 ]), the 
resulting constraint can be written as an LMI constraint. Therefore, Algorithm [1] can be applied to 
solve problem (|5.14[) in the case Di 2 ^ 0. 



5.4 rioo optimization: BMI formulation 

Alternatively, we can also apply Algorithm [1] to solve the optimization with BMI constraints arising 
in Hca optimization of the linear system (II. ip . Let us assume that D 2 i = 0, then this problem is 
reformulated as the following optimization problem with BMI constraints [20] : 



mm 



s.t. 



A T F X + XA F 
BfX 
C F 



C T F ' 



-< 0, X y 0, 7 > 0. 



(5.16) 



Here, as before, A F = A + BFC and C F = Ci + D 12 FC. The bilinear matrix term A T F X + XA F 
at the top-corner of the last constraint can be decomposed as (|5.2|) or (|5.3[) . Therefore, we can 
use these decompositions to transform problem (|5.16p into (I3.5p . After linearization, the resulting 
subproblem is also rewritten as a standard SDP problem by applying Lemma 13.21 We omit this 
specification here. 

To determine a starting point, we perform Phase 1 which is similar to the one carried out in the 
%2-optimization subsection. 
Phase 1. (Determine a starting point 

Step 1. If a (A T + A) < then set F° = 0. Otherwise, go to Step 3. 

Step 2. Solve the following optimization with LMI constraints 



min 7 s.t. 



A T FO X + XA F0 XBi 



BfX 

C f° 



<^f° 



-< 0, X y 0, 7 > 0, 



(5.17) 



where A F o := A + BF°C and C F o := C\ + D 12 F°C. If this problem has a solution 7 and X° 
then terminate Phase 1 and using F° together with 7 , X° as a starting point x° for Phase 2. 
Otherwise, go to Step 3. 



Step 3. Solve the following feasibility problem of LMI constraints: 



Find P y 0, 7 > and K such that: 



PA T +AP+K T B T +BK Bi Pd + K^ 



CiP+D l2 K 



Dn 



^0, 



to obtain K* , 7* and P* . Compute F* := K*(P*)~ 1 C + , where C + is a pseudo-inverse of C, 
and resolve problem (|5.17p with F a :— F* . If problem (|5.17p has a solution X° and 7 then 
terminate Phase 1. Set x° :— (F°, X°,j°). Otherwise, perform Step 4. 



Step 4- Apply the method in Subsection 15.21 to solve the following BMI feasibility problem: 

Find F and P y such that : (A + BFC) T P + P(A + BFC) -< 0. (5.18) 

If this problem has a solution F° then go back to Step 2. Otherwise, declare that no strictly 
feasible point for (|5.16l) is found. □ 
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As in the Tt 2 problem, Phase 1 of the Hoo is also terminated after finitely many iterations. In this 
subsection, we also test this algorithm for several problems in COMPl ib using the same parameters 
and the stopping criterion as in the previous subsection. The computational results are shown in 
Table |3l The numerical results computed by HIFOO and PENBMI are also included in Table [3] 



Table 3: T-Loo synthesis benchmarks on COMPl c ib plants 



Problem 


Other Results, 'Hoc 


Results and Performances 


Name 








n 2 




H1EOO 


PENBMI 




iter 


timeLsJ 


AC1 


5 




3 


2 


3 


0.0000 


_ 


0.0177 


93 


28.050 


A(J2 


5 


3 


3 


5 


3 


0.1115 


- 


0.1140 


99 


32.540 


A(J3 


b 


4 


2 


b 


b 


4.7021 


- 


3.4859 


210 


76.170 


AC4 


4 




1 


2 


•2 


0.9355 


- 


69.9900 


2 


•2.620 


A(J6 


V 


4 


2 


V 


V 


4.1140 


- 


4.1954 


167 


138.570 


ACT 


9 


2 


1 


1 


4 


0.0651 


0.3810 


0.0548 


300 


278.330 


A(J8 


9 


b 


1 


2 


10 


2.0050 


- 


3.0520 


247 


298.070 


AC9 6 


10 


5 


4 


2 


10 


1.0048 


_ 


0.9237 


300 


470.910 


ACll 6 


5 


4 


2 


5 


5 


3.5603 




3.0104 


68 


60.260 


AC12 


4 


4 


3 


1 


3 


0.3160 




2.S02S 


300 


181 870 


A(J15 


4 


a 


2 


6 


4 


15.2074 


427.4106 


15.1995 


105 


36.700 


A(J16 


4 


4 


2 


(j 


4 


15.4060 




14.9881 


186 


H8 82(1 


AC17 


4 


2 


1 


4 


4 


6.6124 




6.6373 


129 


42.400 


HEl b 


4 


1 


2 


2 




0.1540 


1.5258 


0.1807 


300 


143.520 


— in.-i 1 

riE2 


1 

4 


— T — 

2 


— n — 1 

2 


T 1 

4 


T 1 

4 


71 — T7TT5T — 1 

4.4931 




o.784b 


177 


DY.47U 


HE3 


8 


6 


4 


10 


1 


0.8545 


1 .6843 


0.9243 


105 


95.UUU 


HE4" 


8 


6 


4 


12 


8 


23.3448 


- 


22.8713 


252 


325.580 


HE5 6 


8 


2 


4 


4 


3 


8.8952 


- 


37.3906 


300 


430.820 


— IJ L 1 A 1 1 


— a — 1 

4 


— ~> — 1 


— n — 1 

2 


— 2 — 1 

4 


— 2 — 1 

4 


— a aiwa — 1 
0.8975 




— a smi — 
0.8815 


T7T' — 1 

96 


n — 1 ■ v 1 — 

34.4-iO 


REA2" 


4 


2 


2 


4 


4 


1.1881 


- 


1.4188 


300 


118.320 


KEA3 


12 


3 


1 


12 


12 


74.2513 


74.4460 


74.5478 


2 


234.800 


D1S1 


8 


4 


4 


8 


1 


4.1716 




4.1943 


93 


66.130 


D1S2 


3 


2 


2 


3 


3 


1.0548 


1.7423 


1.1546 


54 


17.120 


D1S3 


6 


4 


4 


6 


6 


1.0816 




1.1382 


285 


195.960 


D1S4 


6 


6 


4 


6 


6 


0.7465 




0.7498 


126 


93.220 


TGI 6 


10 


2 


2 


10 


10 


12.8462 




12.9336 


45 


84.380 


AGS 


12 


2 


2 


12 


12 


8.1732 


188.0315 


8.1732 


24 


55.290 


WEC2 


10 


4 


3 


10 


10 


4.2726 


32.9935 


6.6082 


300 


476.010 


WEC3 


10 


4 


3 


10 


10 


4.4497 


200.1467 


6.8402 


300 


425.330 


BD'1'1 


11 


3 


3 


6 


1 


0.2664 




0.8562 


29 


40.910 


MFP 


4 


2 


3 


4 


4 


31.5899 




31.6079 


171 


57.430 


1H 


21 


10 


11 


11 


21 


1.9797 




1.1858 


114 


1206.340 


CSE1 


20 


10 


2 


12 


1 


0.0201 




0.0220 


3 


20.250 


PSM 


7 


3 


2 


5 


2 


0.9202 




0.9227 


87 


67.470 


EB1 


10 


1 


1 


2 


2 


3.1225 


39.9526 


2.0276 


300 


295.420 


EB2 


10 


1 


1 


2 


2 


2.0201 


39.9547 


0.8148 


84 


73.770 


EB3 


10 


1 


1 


2 


2 


2.0575 


3995311.0743 


0.8153 


84 


75.820 


NN1 


3 


2 


1 


3 


3 


13.9782 


14.6882 


18.4813 


300 


144.940 


NN2 


2 


1 


1 


2 


2 


2.2216 




2.2216 


9 


4.060 


NN4 


4 


3 


2 


4 


4 


1.3627 




1.3802 


156 


51.480 


NN8 


3 


2 


2 


3 


3 


2.8871 


78281181.1490 


2.9345 


180 


51.830 


NN9" 


5 


2 


3 


4 


2 


28.9083 




32.1222 


300 


129.920 


NNll" 


16 


5 


3 


3 


3 


0.1037 




0.1566 


9 


366.890 


NN15 


3 


2 


2 


4 


1 


0.1039 




0.1194 


6 


3.740 


INN 16 


8 


4 


4 


4 


8 


0.9557 




0.9656 


48 


37.950 


NN17 


3 


1 


2 


2 


1 


11.2182 




11.2381 


117 


32.160 



Here, the notation is the same as in Table [3J except that denotes the 'H o-norm of the closed- 
loop system for the static output feedback controller. We can see from Table [3] that the optimal 
values reported by Algorithm [T] and HIFOO are almost similar for many problems whereas in general 
PENBMI has difficulties in finding a feasible solution. 



5.5 T-iijT-Loo optimization: BMI formulation 

Motivated from the H2 and Jioo optimization problems, in this subsection we consider the mixed 
H 2 /'Hoo synthesis problem. Let us assume that £>n = 0, D 2 i = and the performance output z is 
divided in two components, z\ and z 2 . Then the linear system (jl.ip becomes: 

x = Ax + B\w + Bu, 
Zl = C?x + Dl 2 u, 

z 2 = C 7 i 2 x + Dgu, { ° y> 

y = Cx. 

The mixed T~L 2 /'H 00 control problem is to find a static output feedback gain F such that, for u = Fy, 
the %2-norm of the closed loop from w to z 2 is minimized, while the "Hoo-norm from w to z\ is less 
than some imposed level 7 [2lJ [27] . 
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This problem leads to the following optimization problem with BMI constraints |27j : 



min trace(Z) 

F,P lt P 2 ,Z 



s.t. 



A T F P 2 



BfPx 
P 2 A F P 2 B l 



(cffcf 



BTPs 



-I 



-< o, 



P\B\ 

-1 2 I 
Pi 



(cpf 
z 



(5.20) 



>- 0, Pi >- 0, P 2 >- 0, 



where A F := A + BFC, Cp 1 := Cf + D^FC and C Z F 2 := Cf + Df%FC. Note that if (7 = J ns> , the 
identity matrix, then this problem becomes a mixed Ti. 2 /Tioa of static state feedback design problem 
considered in [27] ■ In this subsection, we test Algorithm [T] for the static state feedback and output 
feedback cases. 

Case 1. The static state feedback case (C = I nm )- First, we apply the method in [27] to find an 
initial point via solving two optimization problems with LMI constraints. Then, we use the same 
approach as in the previous subsections to transform problem (|5.20p into an optimization problem 
with psd-convex-concave matrix inequality constraints. Finally, Algorithm [1] is implemented to solve 
the resulting problem. For convenience of implementation, we introduce a slack variable 77 and then 
replace the objective function in (|5.16l) by f(x) — rj 2 with an additional constraint trace(.Z') < rj 2 . 

In the first case, we test Algorithm Q] with three problems. The first problem was also considered 
in [H] with 





"-1.40 


-0.49 


-1.93" 




"-0.16 


-1.29~ 




0.25" 


A = 


-1.73 


-1.69 


-1.25 


, B x = 


0.81 


0.96 


, B = 


0.41 




0.99 


2.08 


-2.49 




0.41 


0.65 




0.65 



Cf = [-0.41 0.44 0.68], Cf = [-1.77 0.50 



-0.40 



u \1 



n Z2 



1, and 7 = 2. 



If the tolerance e = 10~ 3 is chosen then Algorithm Q] converges after 17 iterations and reports the 
value rj = 0.7489 with F = [1.9485 0.3990 -0.2119] . This result is similar to the one shown in [2"7] . 
If we regularize the subproblem (|4.1[) with p = 0.5 x 10 -3 and Q = I FF then the number of iterations 
is reduced to 10 iterations. 

The second problem is DIS4 in COMPLfib [H| . In this problem, we set C{ x = Cf and D\\ = £>g 
as in |27) . Algorithm [T] converges after 24 iterations with the same tolerance e = 10~ 3 . It reports 
7] = 1.6925 and 7 = 1.1996 with 



F 



-0.8663 -0.6504 

0.1591 -0.4941 

-0.7017 -0.0785 

-0.0522 -0.5556 



-1.1115 
-0.6322 
0.6121 
-0.5838 



-0.1951 
-0.5409 
-0.8919 
0.4497 



-0.6099 0.2065 ' 

-1.2895 0.2774 

0.2518 -0.2354 

-1.4279 -0.6677 



If we regularize the subproblem (|4.ip with p = 0.5 x 10 and Q = Ip F then the number of iterations 
is 18 iterations. 

The third problem is AC16 in COMPl ib [5D]. In this example we also choose Cf = Cf and 
D\\ = D\\ as in the previous problem. As mentioned in [37], if we choose a starting value 70 = 100, 
then the LMI problem can not be solved by the SDP solvers (e.g., Sedumi, SDPT3) due to numerical 
problems. Thus, we rescale the LMI constraints using the same trick as in [27]. After doing this, 
Algorithm[T]converges after 298 iterations with the same tolerance e = 10~ 3 . The value of 77 reported 
in this case is 77 = 12.3131 and 7 = 20.1433 with 



-1.8533 
4.2672 



0.1737 
-0.9668 



0.6980 
-1.5952 



6.4208 
-2.9240 



The results obtained by Algorithm [T] for solving problems DIS4 and AC16 in this paper confirm the 
results reported in [2"7] , 

Case 2. The static output feedback case. Similarly to the previous subsections, we first propose a 
technique to determine a starting point for Algorithm [1] We described this phase algorithmically as 
follows. 

Phase 1. {Determine a starting point x°). 
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Step 1. If a (A T + A) < then set P° = 0. Otherwise, go to Step 3. 
Step 2. Solve the following linear SDP problem: 



min trace(Z) 



s.t. 



A T p0 P 2 



BJP, 



BfP 1 

P2A F o 
2 



F° 



P2B1 

-I 



■<0, 



P 1 B 1 

-1 2 I 

P2 

L -' po 



^0, 

(c z F l) T 
z 



(5.21) 



y 0, Pi >- 0, P 2 >- 0, 



where A F o = A 



PP°C\ C^o = 



Cf + D^P°C and C* 2 = Cf 



optimal solution Pi,P 2 and Z° then terminate Phase 1. 
point of Algorithm [TJ Otherwise go to Step 3. 
Step 3. Solve the following LMI feasibility problem: 



P* 2 P°C. If this problem has an 
Set x° := (F ,P^P^Z ) for a starting 



Find Q^0,W and Z such that: 
~AQ + QA T +BW+W T B T B x 
B\ -I u 
d+D 12 W O 
AQ + QA T + BW + W T B T 



(Ci 



Pi 

if I) 



O 

-l 2 h 

-< 



dQ + D 12 W 



(CiQ 



- D 12 Wf 
Z 



^0, 



to obtain a solution Q* , W* and Z*. Set P 
C. Solve again problem (|5.21j) with P° := P 
1. Otherwise, perform Step 4. 

5iep ^. Solve the following optimization with BMI constraints 



= W*(Q*) 1 C + , where C + is the pseudo- inverse of 
If problem (|5.21[) has solution then terminate Phase 



max/3 

B,F,P 1 ,P 2 

s.t. Pi >- 0, P 2 >- 0, 

AlP^P.Al + iC^fC^ 



A T F P 2 



± F' 

P 2 A F 



F ) U F 1 7 



iPiPiPfPi ^ -2/3Pi, 



(5.22) 



P 2 B 1 B?P 2 * -2/3P 2 



to obtain an optimal solution P* corresponding to the optimal value /3*. If /3* > then set P° := P* 
and go back to Step 2 to determine Pf, P° and Otherwise, declare that no strictly feasible point 
of problem (|5.20[) is found. □ 

Since at Step 4 of Phase 1, it requires to solve an optimization problem with two BMI constrains. 
This task is usually expensive. In our implementation, we only terminate this step after find a 
strictly feasible point with a feasible gap 0.1. If matrix C is invertible then the matrix P* at Step 
3 is P* = W*(Q*)~ 1 C~ 1 . Hence, we can ignore Step 4 of Phase 1. 

To avoid the numerical problem in Step 3, we can reformulate problem (|5.5[) equivalently to the 
following one: 



Find Q >- 0, W and Z such that: 
~AQ + QA T +BW+W T B T B Y 

Pf -7ha 

Ci + D l2 W O 



(Ci 



-D 12 Wf 

o 



-<o, 



dQ + D 12 W Z 



>~0. 



We test the algorithm described above for several problems in COMPl c ib with the level values 
7 = 4 and 7 = 10. In these examples, we assume that the output signals z\ = z 2 . Thus we have 
C\ x = Ci 2 = Ci and = — D\ 2 . The parameters and the stopping criterion of the algorithm 
are chosen as in Subsection 15.31 The computational results are reported in Table [4] with 7 = 4 and 
7 = 10. Here, 'H 2 /'H ao are the *H 2 and Hoo norms of the closed- loop systems for the static output 
feedback controller, respectively. With 7 = 10, the computational results show that Algorithm [1] 
satisfies the condition ||Poo(s)||oo < 7 = 10 for all the test problems. While, with 7 = 4, there are 
5 problems reported infeasible, which are denoted by The 'Hoo-constraint of three problems: 
AC3, AC 11 and NN8 is active with respect to 7 = 4. 
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Table 4: T^/'Hoo synthesis benchmarks on COMPl G ib plants 



Problem 


Results and Pcrlormanecs (7 — 4) 


Results and Pcrlormanecs (7 — 10) 


Name 




iter 


time[s] 




iter 


timeCs] 


AC1 


0.0585 / 0.0990 


3 


4.22 


0.0585 / 0.0990 


3 


4.27 


AC 2 


0.1067 / 0.1723 


6 


7.31 


0.1070 / 0.1727 


3 


7.15 


AC 3 


5.2770 / 3.9999 


51 


281.53 


4.5713 / 5.1298 


18 


19.18 


AC 6 


- / - 


- 


- 


4.0297 / 4.8753 


283 


330.64 


AC 7 


0.0415 / 0.0961 


1 


3.39 


0.0420 / 0.1286 


2 


3.91 


AC 8 


1.2784 / 2.2288 


43 


60.78 


1.3020 / 2.5719 


23 


31.59 


AC11 


4.0704 / 4.0000 


76 


175.75 


4.0021 / 5.1949 


117 


122.86 


AC12 


0.0924 / 0.3486 


18 


73.46 


1.4454 / 1.6444 


300 


234.13 


AC17 


- 1 - 


- 


- 


4.1228 / 6.6472 


2 


11.620 


HE1 


0.1123 / 0.2257 


2 


131.18 


0.0973 / 0.2080 


1 


30.97 


HE2 


- / - 


- 


- 


4.7302 / 9.8931 


75 


55.48 


REA1 


1.8214 / 1.4740 


30 


25.64 


1.8213 / 1.4730 


30 


26.65 


REA2 


3.5014 / 3.5180 


42 


22.09 


3.5015 / 3.5209 


45 


23.26 


DIS1 


- / - 






2.8505 / 4.7904 


15 


30.51 


DIS2 


1.5079 / 1.8500 


18 


7.92 


1.5079 / 1.8520 


21 


7.92 


DIS3 


2.0577 / 1.7934 


27 


25.03 


2.0577 / 1.7766 


30 


24.54 


D1S4 


1.6926 / 1.1952 


21 


18.62 


1.6926 / 1.2009 


24 


21.55 


AGS 








7.0332 / 8.2035 


8 


196.73 


PSM 


1.5115 / 0.9248 


177 


160.41 


1.5115 / 0.9248 


180 


167.31 


EB2 


0.7765 / 1.0828 


7 


9.70 


0.7768 / 1.0867 


10 


13.16 


EB3 


0.8406 / 0.9249 


1 


3.21 


0.8383 / 0.9418 


1 


2.93 


EB4 


1.0147 / 1.0707 


6 


59.55 


0.9981 / 1.2146 


12 


111.26 


NN2 


1.5651 / 2.4834 


12 


5.37 


1.5651 / 2.4876 


12 


5.49 


NN4 


1.8778 / 2.0501 


202 


154.49 


1.8779 / 2.0519 


213 


161.00 


NN8 


2.3609 / 3.9999 


21 


15.71 


2.3376 / 4.6514 


15 


6.57 


NN15 


0.0820 / 0.1010 


42 


18.75 


0.0771/0.1012 


24 


10.47 


NN16 


0.3187 / 0.9574 


90 


96.44 


0.3319 / 0.9572 


258 


303.87 



6 Concluding remarks 

We have proposed a new algorithm for solving many classes of optimization problems involving BMI 
constraints arising in static feedback controller design. The convergence of the algorithm has been 
proved under standard assumptions. Then, we have applied our method to design static feedback 
controllers for various problems in robust control design. The algorithm is easy to implement using 
the current SDP software tools. The numerical results are also reported for the benchmark collec- 
tion in COMPl e ib. Note, however, that our method depends crucially on the psd-convex-concave 
decomposition of the BMI constraints. In practice, it is important to look at the specific struc- 
ture of the problems and find an appropriate psd-convex-concave decomposition for Algorithm [TJ 
The method proposed can be extended to general nonlinear semidefinite programming, where the 
psd-convex-concave decomposition of the nonconvex mappings are available. From a control design 
point of view, the application to more general reduced order controller synthesis problems and the 
extension towards linear parameter varying or time-varying systems are future research directions. 
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Proof of Lemma 14.21 

For any matrices A, B £ S+, we have A o B > 0. From Step 1 of Algorithm [TJ we have x k+1 is the 
solution of the convex subproblem (|4.1j) and A fc+1 is the corresponding multiplier, under Assumption 
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x k+1 ) 



[31 they must satisfy the following generalized Kuhn- Tucker condition: 

G Wf(x k ^)+p k QlQ k (x^ 1 -x k ) + {Y:UiD[G i (x)-H i {x k ) 

-DHitfKx - x k )]\ xk+ ,y A k+1 + N n (x k+1 ), ( 1} 

Gi{x k+1 ) - H t {x k ) - DH{x k ){x k+1 - x k ) < 0, A* h 0, 
[G i (x k+1 )-H i (x k )-DH{x k )(x k+1 -x k )] oA k+1 = 0. 

Noting that D [Gi(x] ~ H l {x k ) - DHi(x k )(y - x k )] | x=x *+i = DG l (x k+1 ) - DH l (x k ) for i = 1, . . .1, 
it follows from the first line of Q| and the convexity of / that 

/(y)-/(^ +1 )+|E^(^ +1 )-^(^)]*A" +1 | (y-* k+1 ) 

> l^fix^+J^DG^+^-DH^WA^ {y 

+ f lly-^ll^f h-x^g+Pkiy-x^fQlQ^-x^), Vj/GO. (.2) 
On the other hand, we have 

{[DG l {x k+1 ) - DH,(a; fc )]*A* +1 } T (y - x fc+1 ) 

= A* +1 o [DGi(a; fc+1 )(y - x k+1 ) - DH t (x k )(y - a: fe+1 )]. (.3) 

Since Gi and 7?i are psd-convex, applying Lemma 12. II we have 

Gi(x k )-Gi(x k+1 )>DGi (x k+1 ) (x k ~x k+1 ), 
and H l (x k+1 )- H l (x k )hDH l (x k )(x k+1 -x k ), i = l,...,l. 

Summing up these inequalities we obtain 

G l (x k ) - Hi (x k ) - [d (x k+1 ) - Hi (x k+1 )] h [Dd (x k+1 ) (x k - x k+ 1 ) - DH t {x k ) (x k - x k+ 1 )] . 

Using the fact that A k+1 >z 0, this inequality implies that 

A k+1 o{Gi{x k )- H^) - [G t (x k+1 ) - H^ 1 )]} 

h A k+1 o [DG t (x k+1 )(x k - x k+1 ) - DH t (x k )(x k - x k+1 )}. (.4) 

Substituting y = x k into ([21) and then combining the consequence, Q5]), ([31) and the last line of (p"j) 
to obtain 

z 

f(x k )-f(x M1 )+Y l A? 1 °[Gi{x'>)-H i (x k )] > ^\\x k+1 -x k \\l + Pk \\Q k {x k+1 -x k )\\l (.5) 

i=l 

Now, since x k is the solution of the convex subproblem (|4.1[) linearized at a:' 0-1 . One has Gi{x k ) — 
Hi(x k ) < 0. Moreover, since A k+1 t 0, we have A^' +1 o [Gi(x k ) - Hi(x k )] < 0. Substituting this 
inequality into (0, we obtain 

f(x k ) - f(x k+1 ) > ^\\x k - x k+l \\ 2 2 + p k \\Q k (x k+1 - x k )\\l 

This inequality is indeed (|4.3|) which proves the item a). If there exists at least one io G {1, . . . , 1} 
such that G io (x k ) - H l0 (x k ) -< and A^ 1 y then A^ 1 o [G ia {x k ) - H ia {x k )] < 0. Substituting 
this inequality into Q2) we conclude that f(x k+1 ) < f{x k ) which proves item b). The last statement 
c) follows directly from the inequality (|4.3p . □ 
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